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Abstract 

Reconstruction of the curvatures of radio wavefronts of air showers initiated by 
ultra high energy cosmic rays is discussed based on minimization algorithms 
commonly used. We emphasize the importance of the convergence process in- 
duced by the minimization of a non-linear least squares function that affects the 
results in terms of degeneration of the solutions and bias. We derive a simple 
method to obtain a satisfactory estimate of the location of the main point of 
emission source, which mitigates the problems previously encountered. 

Keywords: UHECR, radio-detection, antennas, non-convex analysis, 
optimization, ill-posed problem. 



1. Introduction 

The determination of the nature of the ultra-high energy cosmic rays (UHECR) 
is an old fundamental problem in cosmic rays studies. Numerous are the difficul- 
ties. New promising approaches could emerge from the use of the radio-detection 
method which exploits, through antennas, the radio signal that accompanies the 
development of the extensive air shower (EAS). Several experimental prototypes 
like CODALEMA [1 in France and LOPES [3 in Germany shown the feasibility 
and the potential of the method to reconstruct EAS parameters, as the arrival 
direction, the impact location at ground, the lateral distribution function of 
the electric field, or the primary particle energy jH [H [7j |8]. However, the 
temporal radio wavefront characteristics remain still poorly determined |51 ITU], 
although its knowledge could be consider as one of the first steps in retrieving 
information about the EAS itself. The importance of this information resides 
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in its potential sensitivity to the nature of the primary particle, especially be- 
cause the existence of a curvate radio wavefront (a spherical wavefront) could 
provide the location of the main point of the emission source, and possibly an 
estimation of Xmax, event by event. Indeed, the arrival timing being defined by 
the maximum amplitude of the radio signal, it is more likely linked to a limited 
portion of the longitudinal development of the shower (and so especially at the 
point of maximum) . 

Moreover, the migration of present small scale radio-prototypes to large scale 
experiments spread over surfaces of several tens of 1000 km 2 using self-triggered 
antennas, is challenging. This technique is subjected to delicate limitations 
in regard to UHECR recognition, due to noises induced by human activities 
(high voltage power lines, electric transformers, cars, trains and planes) or by 
stormy weather conditions (lightning) . Figure [I] shows a typical reconstruction 
of sources obtained with the CODALEMA experiment [12], by invoking a spher- 
ical wave minimization. Such patterns are also commonly observed in others 
radio experiments [HI [IB] . In most of the cases, one of the striking results is that 
these emission sources are reconstructed with great inaccuracy, although they 
are fixed and although the number of measured events is high. By extension, a 
cosmic event being a single realization of the detected observables (arrival time 
and peak amplitude on each antenna), interpretation of such methods of recon- 
struction for the identification of a point source can become even more delicate, 
even using statistical approaches. 
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Figure 1: Typical result of reconstruction of two entropic emitters at ground, observed with the 
stand-alone stations of CODALEMA, through standard minimization algorithms. Despite the 
spreading of the reconstructed positions, these two transmitters are, in reality, two stationary 
point sources. 



The commonly used technique relies on the minimization of an objective 
function which depends on the assumed shape of the wavefront, using the arrival 
times and locations of the antennas. The aim of this paper is to highlight 
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that the minimization of such an objective function, incorporating a spherical 
wave front, can be an ill-posed problem. We will show that it originates from 
strong dependencies of the convergence of the minimization algorithms with 
initial parameters, from the existence of degenerations of the solutions (half 
lines) which can trap most of the common algorithms, and from the existence of 
offsets (bias) in the reconstructed positions. Finally, by avoiding more complex 
estimates based on advanced statistical theories, we got to deduce a simple 
method to obtain a significant estimate of the source location. We compared 
the exact results with our numerical reconstructions performed on a test array. 



2. Reconstruction with common algorithms 

The performances of different algorithms has been tested using the simplest 
test array of antennas. Within the constraints imposed by the number of free 
parameters used for reconstruction, we choose an array of 5 antennas for which 
the antennas positions rf = (xi,yi,Zi) are fixed (see Fig. |2| (this corresponds 
to a multiplicity of antennas similar to that sought at the detection threshold 
in current setups). 
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Figure 2: Scheme of the antenna array used for the simulations. The antenna location is took 
from a uniform distribution of 1 m width. 



A source S with a spatial position r s = (x s ,y s ,z s ) is set at the desired 
value. Assuming t s the unknown instant of the wave emission from S, c the 
wave velocity in the medium considered constant during the propagation, and 
assuming that the emitted wave is spherical, the reception time fj on each 
antenna i G {1, . . . , N} can written: 



\J (xi - x s f + (yi - y s ) 2 + (z~i - z s y 



+ G{0,a t ) 
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where G(0, at) is the Gaussian probability density function centered to t = 
and of standard deviation at- This latter parameter stand for the the global 
time resolution, which depends as well on technological specifications of the 
apparatus than on analysis methods. 

The theoretical predictions are compared to the reconstructions given by 
the different algorithms. The latter are setup in two steps. First, a planar 
adjustment is made, in order to pres-tress the region of the zenith angle and 
azimuth angle ef> of the source arrival direction. It specifies a target region in 
this subset of the phase space, reducing the computing time of the search of the 
minimum of the objective function of the spherical emission. Reconstruction 
of the source location is achieved, choosing an objective-function that measures 
the agreement between the data and the model of the form, by calculating the 
difference between data and a theoretical model (in frequentist statistics, the 
objective-function is conventionally arranged so that small values represent close 
agreement): 

1 N 2 

/(^ ! *:) = 2E[ii^-^n 2 -^-**) 2 ] c 1 ) 

The partial terms ||rt — rt\\ 2 — (t* — t*) 2 represents the difference between 
the square of the radius calculated using coordinates and the square of the 
radius calculated using wave propagation time for each of the N antennas. The 
functional / can be interpreted as the sum of squared errors. Intuitively the 
source positions rt at the instant t s is one that minimizes this error. 

In the context of this paper, we did not use genetic algorithms or multivari- 
ate analysis methods but we focused on three minimization algorithms, used 
extensively in statistical data analysis software of high energy physics [T"5l I16j : 
Simplex, Line-Search and Levenberg-Marquardt (see table [2J. They can be 
found in many scientific libraries as the Optimization Toolbox in Matlab, the 
MPFIT in IDL and the library Minuit in Root that uses 2 algorithms Migrad and 
Simplex which are based respectively on a variable-metric linear search method 
with calculation of the objective function first derivative and a simple search 
method. For the present study, we have used with their default parameters. 

We tested three time resolutions with times values took within 3 at ■ 

• at = Oris plays the role of the perfect theoretical detection and serves as 
reference; 

• at = ins reflects the optimum performances expected in the current state 
of the art; 

• at = 10 ns stands for the timing resolution estimate of an experiment like 
CODALEMA [2]. 

For every source distance and temporal resolution, one million events were gen- 
erated. Antenna location was taken in a uniform distribution of 1 m width. A 
blind search was simulated using uniform distribution of the initial r s values 
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Table 1: Summary of the different algorithms and methods used to minimize the objective- 
function. The second row indicates framework functions corresponding to each algorithm; 
third recalls the framework names. The key information used for optimization are recalled 
down, noting that a differentiable optimization algorithm (ie. non-probabilistic and non- 
heuristic) consists of building a sequence of points in the phase space as follows: Xk+i = 
x k an d that it is ranked based on its calculation method of and parameters (see 

El). 
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from 0.1 km to 20 km. Typical results obtained with our simulations are pre- 
sented in Figures [3] and |4j The summary of the reconstructed parameters is 
given in table [2] 



^=1 km, HT, « LVM algorithm ^ k ™. »° W ™ algorithm 




Figure 3: Results of the reconstruction of a source with a radius of curvature equal to 1 and 
10 km with the LVM algorithm. For Rtrue = 1 km, the effect of the blind search leads to non- 
convergence of the LVM algorithm, when initialization values are greater than Rt rU e = 1 km. 



R Tm =1 km, 8=90", M°, Simplex algorithm 11=10 km, 6=90°, (rf, Simplex algorithm 




Figure 4: Results of the reconstruction of a source with a radius of curvature equal to 1 and 
10 km with the Simplex algorithm. 

Whatever the simulations samples (versus any source distances, arrival direc- 
tions, time resolutions), (also with several detector configurations) and the three 
minimization algorithms, large spreads were generally observed for the source 
locations reconstructed. This suggests that the objective-function presents local 
minima. Moreover, the results depend strongly on initial conditions. All these 
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Reconstruction Results 


a t (ns) 


Rtrue{m) 


Algorithms 


Rmean (^) 


Rmodeim) 


a R (m) 





1000 


Levenberg-Marquardt 


(10071) 1002 


(1081) 1081 


(5763) 102 


Simplex 


1198 


1133 


1477 


3000 


Levenberg-Marquardt 


(9960) 3082 


(2998) 2998 


(5781) 302 


Simplex 


3134 


3272 


3437 


10000 


Levenberg-Marquardt 


9999 


9997 


56 


Simplex 


10466 


9929 


5817 


3 


1000 


Levenberg-Marquardt 


(10071) 1003 


(934) 934 


(5763) 108 


Simplex 


1199 


168 


1486 


3000 


Levenberg-Marquardt 


(9954) 3068 


(2874) 2874 


(5792) 495 


Simplex 


3132 


3010 


3485 


10000 


Levenberg-Marquardt 


7174 


6877 


3021 


Simplex 


8194 


6479 


6154 


10 


1000 


Levenberg-Marquardt 


(10068) 985 


(964) 964 


(5767) 175 


Simplex 


1189 


199 


1507 


3000 


Levenberg-Marquardt 


(9703) 2238 


(1620) 1620 


(6125) 877 


Simplex 


2760 


2070 


3703 


10000 


Levenberg-Marquardt 


2770 


667 


2305 


Simplex 


3675 


934 


4048 



Table 2: Summary of parameters reconstructed with different algorithms for several distances 
of source and several timing resolutions. On the Levenberg-Marquardt, the results in paren- 
theses are those taking into account the flat portion of the resulting distribution (see Fig. 
[3]). They are typical of initialization values which are starting too far from the actual source 
distance. The Line-Search method was ultimately rejected for this quantitative study, because 
results too dependent on the starting algorithm fixing the initial conditions. 
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Figure 5: Condition numbers obtained using the formula Cond(Q) = |]Q||.||Q — 1 ||with Q the 
Hessian matrix (see next section) as a function of the source distance and for different timing 
resolutions. The large values of conditioning suggest that we face an ill-posed problem. 



phenomena may indicate that we are facing an ill-posed problem. Indeed, con- 
dition number calculations [H] (see Fig. [2j, which measures the sensitivity of 
the solution to errors in the data (as the distance of the source, the timing res- 
olution, etc.), indicate large values (> 10 4 ), when a well-posed problem should 
induce values close to 1. 

To understand the observed source reconstruction patterns, we have under- 
taken to study the main features of this objective-function. 



3. Study of the objective-function for the spherical emission 

To estimate the source position X s — (f s ,t s ) using the sequence of arrival 
ti, the natural method is to formulate an unconstrained optimization problem 
of type a non- linear least square |17| . starting from eq. 1. which can rewrit^] 
(see the notations listed in table [3]): 

1 N 2 1 N 

f (x s ) = -j2[\\T$-rt\\l- (t; - 1*) 2 ] =-j2 ft{x s ) (2) 

i=l i=l 

Several properties of the objective-function / were studied: the coercive 
property to indicate the existence of at least one minima, the non-convexity 
to indicate the existence of several local minima, and the jacobian to locate 



1 In practice of the minimization, it is usual to take into account errors on the measured 
parameters by putting them in the objective function denominator. In our theoretical study, it 
is assumed that the arrival times errors are the same for all the antennas (at = constant V i). 
The present studied functional is generic and does not include errors, but as will see later, 
introduction of a multiplicative constant doesn't change the results of our study. 



3 . 1 Con vexi ty prop erty 



r~l, fl: position of the source, position of the i th antenna 
t s , ti'. emission time of the signal, signal arrival time at the i th antenna 
t*, t*: reduced time variables (ie. t* — c.t) 
a\: time resolution on the i th antenna 
X s , Xi spacio-temporal position of the source, of the i th antenna 
V/, V 2 /: first and second derivative of the objective function / 



M = I4 — 2E 44 = 







-1 
, L 



: second order tensor related to the Minkowski metric 



quadratic and linear form 
< .|. >: inner product 
transpose of a vector or a matrix X 



Table 3: List of notations 



the critical points. (Bias study, which corresponds to a systematic shift of the 
estimator, is postponed to another contribution). In mathematical terms, this 
analysis amounts to: 

• Estimate the limits of / to make evidence of critical points; obviously, the 
objective function / is positive, regular and coercive. Indeed, / tends to 
+00 when \\X\\ — > ±00, because it is a polynomial and contains positive 
square terms. So, / admits at least a minimum. 

• Verify the second optimality condition: the convexity property of a func- 
tion on a domain for a sufficiently regular function is equivalent to positive- 
definiteness character of its Hessian matrix. 

• Solve the first optimality condition: Vf(X s ) = (jacobian) to find the 
critical points. 

3.1. Convexity property 

Using fi{X s ) = (X s — Xi) T .M.(X S — Xi) where M designates the Minkowski 
matrix and given V fi(X s ) = 2.M(X S —Xi), the / gradient function can written 
(see appendix 1): 

\vf{X s ) = (]T fi(X s ))M.X s - M.(£, fi(Xs)Xi) 
The Hessian matrix, which is the / second derivative can written: 

v 2 ./ (x s ) = ]T vfi(x s ).vfl + £ h^ 2 h 

that becomes, replacing V/j by its expression: 

V 2 f(X s ) = (J2MX s )).M+2M.[N.X s .Xj+J2x i X^-X s (J2x i ) T -(J2x t )Xj}.M 



3.2 Critical points 
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Using a Taylor series expansion to order 2 (see appendix 1), an expanded form of 
the Hessian matrix, equivalent to the previous formula of the / second derivative, 
is: 

' & Ki + 2E,(' S -',) ! >E((««-«i)(»<-V() >E((»s-»i)(«.->i) »£i (»» -■»<)(*<- 'j) 

1 * Ei Ki + 2£i („, - Vi ) 2 2^ (v, - Vi) («. - «<) 2E; (v« - Hi) (t* - tjj 

2 * * Ei K i + 2 E; («« - *i) 2 2£i(>»-«i)('*-<y 

- Ei K< + 2Ei (** - *s) 2 
(3) 

This latter allowed us to study the convexity of / (see appendix 1). Indeed, 
because its mathematical form is not appropriate for a direct use of the convexity 
definition, we have preferred to use the property of semi-positive-definiteness of 
the Hessian matrix. Our calculus lead to the conclusion that: 

• Using the criterion of Sylvester [TS] and the analysis of the principal minors 
of the Hessian matrix , we find that / is not convex on small domains, 
and thus is likely to exhibit several local minima, according to X s and 
Xj. It is these minima, which induce convergence problems to the correct 
solution for the common minimization algorithms. 

3. 2. Critical points 

The study of the first optimality condition (Jacobian = 0) gives the following 
system Vf(X s ) = and allows finding the critical points and their phase-space 
distributions. Taking into account the following expression: 

\Vf{X s ) = (E fi(X s ))M.X s - M.(E fi(X„)Xi) 

we get the relation: 

Y _ fi( x *)_ x / 4 \ 

This formula looks like the traditional relationship of a barycenter. Thus, we 
interpret it in terms of the antennas positions barycenter and its weights. The 
weight function fi expressing the space-time distance error between the posi- 
tion exact and calculated, the predominant direction will be the one presenting 
the greatest error between its exact and calculated position. The antennas of 
greatest weight will be those the closest to the source. 

In practice (see appendix 2) , because the analytical development of this op- 
timality condition in a three-dimensional formulation is not practical, especially 
considering the nonlinear terms, we chose to study particular cases. We consid- 
ered the case of a linear antennas array (ID) for which the optimality condition 
is easier to express with an emission source located in the same plane. This 
approach allows us to understand the origin of the observed degeneration which 
appears from the wave equation invariance by translation and by time reversal 
(known reversibility of the wave equation in theory of partial differential equa- 
tions) and provides us a intuition of the overall solution. It also enlightens the 
importance of the position of the actual source relative to the antennas array 



3.3 Convex hull 
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(the latter point is linked to the convex hull of the antenna array and is the 
object of the next section). Our study led to the following interpretations: 

• The iso-barycenter of the antenna array (of the lit antennas for a given 
event) plays an important role in explaining the observed numerical degen- 
eration. The nature of the critical points set determines the convergence 
of algorithms and therefore the reconstruction result. 

• There are strong indications, in agreement with the experimental results 
and our calculations (for ID geometry), that the critical points are dis- 
tributed on a line connecting the barycenter of the lit antennas and the 
actual source location. We used this observation to construct an alterna- 
tive method of locating the source (section 4). 

• According to the source position relative to the antenna array, the recon- 
struction can lead to an ill-posed or well-posed problem, in the sense of J. 
Hadamard. 

3.3. Convex hull 

In the previous section we pointed that to face a well-posed problem (no 
degeneration in solution set), it was necessary to add constraints reflecting the 
propagation law in the medium, the causality constraints, and a condition link- 
ing the source location and the antenna array, the latter inducing the concept 
of convex hull of the array of antennas. From appendix 2, we also saw that 
analytically the critical points evidence could become very complex from the 
mathematical point of view. Therefore, we chose again an intuitive approach to 
characterize the convex hull, by exploring mathematically the case of a linear 
array with an emission source located in the same plane. This is the subject of 
the appendix 3. 

The results extend to a 2D antenna array, illuminated by a source located 
anywhere at ground, arguing that it is possible to separate the array into sub- 
arrays arranged linearly. The superposition of all the convex segments of the 
sub-arrays leads then to conceptualize a final convex surface, built by all the 
peripheral antennas illuminated (see Figure [6]) . 

The generalization of these results to real practical experience (with a source 
located anywhere in the sky) was guided by our experimental observations (per- 
formed through minimization algorithms) that provide a first idea of what hap- 
pens. For this, we chose to directly calculate numerically the objective function 
for both general topologies: a source inside the antenna array (ie. and at ground 
level) and an external source to the antenna array (in the sky ). As can be de- 
duced from the results (see Figs. [7] and [8]), for a surface antenna array, the 
convex hull is the surface defined by the antennas illuminated. (An extrapo- 
lation of reasoning to a 3D network (such as Ice Cube, ANTARES,...) should 
lead, this time, to the convex volume of the setup). 

Our results suggest the following interpretations: 



3.3 Convex hull 
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Translation along the barycentrie direction 




Antenna 5 



Figure 6: Scheme of the reconstruction problem of spherical waves for our testing array of 
antennas (2D), with a source located at ground. For this configuration, the convex hull 
becomes the surface depicted in red. The result is the same for a source in the sky. 




Figure 7: Plots of the objective- function versus R and versus the phase space (R, t), in the 
case of our testing array (2D), for a source on the ground and located inside the convex surface 
of the antenna array. This configuration leads to a single solution. In this case the problem 
is well-posed. 
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Figure 8: Plots of the objective-function versus R and versus the phase space (R, t), in the 
case of our testing array (2D) for a source outside the convex hull. This configuration leads 
to multiple local minima. All minima are located on the line joining the antenna barycenter 
to the true source. In this case the problem is ill-posed. 

• If the source is in the convex hull of the detector, the solution is unique. In 
contrast, the location of the source outside the convex hull of the detector, 
causes degeneration of solutions (multiple local mininiunis) regarding to 
the constrained optimization problem. The source position, outside or 
inside the array, affects the convergence of reconstruction algorithms. 

4. Proposed method of reconstruction 

Taking into account of the previous results, in order to avoid the trap of 
the local minima with common algorithms, we chose to compute directly the 
values of the objective function on a grid, using a subset of phase space in 
the vicinity of the solution a priori, and assuming that the minimum of the 
objective function correspond to the best estimate of the position of the source 
of emission. The input parameters are set from the planar fit, which provides 
both windows in 8 and </>. By cons, this method, being based on the search of 
the minimum of the objective function using a grid calculation, the choice of 
the metric becomes crucial. On the zenith and azimuth angles, the metric can 
be adequately inferred from the value of the angular resolution obtained by the 
current detectors, or 0.1 for 9 and <fi. Looking at the space metrics, the latter 
can be inferred by considering the quantity (c 2 <7 2 ) 2 , ie. around one meter. The 
direction-priori is given by the planar fit, while the quantity r s is left free in 
the range 0.1 — 20 km (the upper bound being determined by the value of the 
curvature exploitable, given the temporal resolutions currently available). 

A typical result obtained with our method is presented in Fig. [9] and a 
summary of the reconstructed parameters is given in table [4] which have to be 
compared to those presented in the table [2j 
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Figure 9: Histogram of the reconstruction of a source located at 10 km from the detector 
array, and using the grid method associated to the search of the global minimum for each 
event. 
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Reconstruction results 
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Table 4: Summary of parameters reconstructed with our method for the same input values as 
given in table 3. 
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5. Conclusion 

Experimental results indicated that the common methods of minimization 
of spherical wavefronts could induce a mis-localisation of the emission sources. 
In the current form of our objective function, a first elementary mathematical 
study indicates that the source localization method may lead to ill-posed prob- 
lems, according to the actual source position. To overcome this difficulty, we 
developed a simple method, based on grid calculation of the objective function. 
This approach appears to provide, at worst, an estimate as good as for the com- 
mon algorithms for locating the main point of the emission source, keeping in 
mind that this method is not optimal in the sense of optimization theories. How- 
ever, further developments are without any doubt still necessary, maybe based 
on advanced statistical theories, like for instance by adding further information 
(as the signal amplitude or the functional of the radio lateral distribution) . This 
could be achieved by trying a generalized objective function which includes these 
parameters. In addition, the interactions with other disciplines which face this 
problem could also provide tracks of work (especially regarding earth sciences 
which focus on technics of petroleum prospecting) . 

6. Appendix 1 

6.1. Symbolic calculus 

Keeping the same notation as in table [3j the objective function can be writ- 
ten: 

1 N 
/(*«)= 2 

i=l 

with fi (X s ) = {X s - X>f ■ M ■ (X s -Xi) = \\¥t- Ttf - (t* s - t*) 2 . 

The formula V/ (X s ) = Y2 fi (Xs) ' ^ fi (Xs) is derived from the formula of 
a product derivation. Using the bi-linearity of the inner product, we show that 
V/, (X s ) — 2M ■ (X s — Xi). By injecting this formula into the formula of V/ , 
we obtain the following formula: 

Vf(X s ) = ^2fi(X s )-Vfi(X s ) 

= Y,fi{X s )-2M-{X s -X i ) 

It then leads to the following form: 

iv/ (X s ) = (J^ fi (Xs)) M ■ X s - M ■ (J2 fi (X s ) X,) 

With the same method, the second derivative matrix (Hessian matrix) is given 
by the following formula : 

V 2 / (x s ) - v /> (Xs) ■ v/, (x s f +J2fi (x s ) ■ v 2 /, (x s ) 



6.2 Taylor expansion and explicit calculus 



16 



By injecting in the previous formula the following formula of the second 
derivatives V 2 /i (X s ) = 2M and by using the relation (AB) = B T A T , we get 
the following formula: 



2M 



V 2 / (X s ) = J2 V/i (X.) ■ V/, (x s f + (*.) • V 2 /, (X.) 
= J2 2M ■ (X s - Xi) • (2M • (Xs - X^f + fi 
=4M ■ [J2 (X s Xj - X s Xj ~ XiXj + X t X?)] ■ M + 2 (j^ fi (X)) ■ M 

NXsXj + ]r x l xj - x s (J2 x ) T - (E X ) X 



=4M • 



■ M 



Both relationships correspond to the end-calculus forms given in the 3 rd section. 
These forms are easy to handle for symbolic calculus but not convenient for 
explicit calculation used for studying the convexity. 



6.2. Taylor expansion and explicit calculus 

An explicit form for the objective function first and second differential can 
be obtained using a Taylor expansion. Indeed, the function/ is an element of 
C°° (M 4 ,IR) J 2 ] and is therefore differentiable in the sense of Frechet. 



Let X s — (rt, t* ) ± be a fixed vector of R 4 and e" 



— 



7i> 



another vector 



of] 



In order to simplify the calculus, we use the following notations : 
Ki = || rt — rj |L — (t* s — t*) 2 a constant term when setting the vector X s ; 



Li (-?) = (Tt-^\it)-(t* s -n) 

and q(m*) 



t* the linear form: 







)= 


t 



t* 2 the quadratic form. 



The Taylor expansion leads to: 
f(X. 



^E(||^+^-^|| 2 -^+**-**) : 
\T,((^+^-^\^+t-Tt)-(t:+f-t;fy 

i 



\\t + ' + 2 (ft - n 1 1) - (t; - t-f - c 2 - 2? k - (•)) 



Using the multinomial expansion, the function / can then be approximated 
by the second-order Taylor expansion following: 

/ (rt + 1, t* + «*) « K *+ 2 E K i- L * **) + 2 E L2 i 0* > t *)+[Yl Ki ) Q p* > r 

i i i \ i J 



2 The function is also an element of the algebra R [Xi, . . . , X4] 



6.3 Convexity property 
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We identify from this formula: 
the constant term ; 

i 

the linear term which is V/ {X S ) T ■ it - 

f first differential in (rj,t*) ); 

and the quadratic form at the point X s 



t* - 1: 



l T 



■ ~f (the 



-Q (X s ,Xi) = 



Ei & (« 



2Ei ( = s - *<) (»i ~ «») 2E, (<»« - *i) (* 
Ei + 2 2 ; (y s - Vi ) 2 2JZi (Vs - Vi) {* 
Ei Ki + 2Ei (is 



-*♦) 2 Ei(ys-»i)( 

- Ei A"i + 2£i (** 



which is the / Hessian matrix in (r\, , <*), or the second differential of / also 
denoted V 2 / (X s , Xi). The use of * indicates that the coefficients above and 
below the diagonal are equal (Schwarz Lemma). The quadratic form represented 
by this matrix gives us the local second-order properties for the function /. To 
show that a critical point is a local minimum, it will suffice to verify that the 
Hessian matrix is definite positive in the vicinity of this point. 

6.3. Convexity property 

The convex analysis occupies a capital place in the problems of minimization. 
Indeed, an important theorem yet intuitive stated that if a convex function has 
a local minimum, it is automatically global. We will shows that the function / 
is not convex in M 4 , i.e. that the Hessian matrix in non-positive define. 

Let V 2 / (X) the Hessian matrix, and let's suppose d a vector, since the func- 
tion / is twice differentiable, using the Sylvester's criterion [18] to characterize 
the convexity of / , we can write the following equivalence: 

/ is convex Hessian is positive semi-definite All Hessian principal minors are just nonnegative 

/ is convex o Vd, MX, d T ■ V 2 / (X)-d^O 

So if we can find an element X and d such as d T ■ V 2 / (X) ■ d < 0, / will be 
non-positive definite. For this, it is sufficient to find a single negative principal 
minor to demonstrate the Hessian matrix is non-positive definite. The objective 
function f will present then several local minimums and will be thus locally non- 
convex. 

So let Q the explicit expression of the Hessian and let us choose d T = (0001) 
then: 

/°\ 

d T ■ V 2 f{X) ■ d = (0 1) • Q(X S , X^ ■ 






V i J 



i i 

which is represent the principal minor of order 4 of the Hessian. 
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For a family of fixed positions antennas and for a signal source with coor- 
dinates X s such as y s — z s = t* = 0, the negativity condition of the principal 
minor of order 4 can then written: 

i i 

Now the left term tends to infinity when the source tends to infinitjj^] It is 
written in terms of limits, 

lim (x s - x % ) 2 = +oo VA > 0, 3rj > 0/\x s \ > rj => (x a - xA 2 > A 

\x s |— f+oo ^ — ' ^ — ' 
i i 

Taking a value {~Ui — z f + 3if*) of the constant A, it exist a real 77 and 

therefore a a: s such that (x s — xi) 2 > ^ (— y| — + 3i* 2 ). We deduce that 

the function is not convex in the vicinity of this point. It suffices to take d T = 
(0 1) and x s =r/ + l. 



7. Appendix 2 

7.1. Degeneration line for a linear antenna array 

According to experimental data analysis and to our simulations (see Fig. [I] 
and [8| , the results of the common minimization algorithms appear to fall on a 
half-line in the phase space (x, y, z) which we shall call the degeneration line, 
which is linked to the existence of local minima. We present the mathematical 
development in the case of a linear array using an analysis-synthesis method. 
Then we try to generalize results to the higher dimension cases. 

Let suppose X s = (x s , t* s ) a critical point of /, ie. V f(X s ) = 0, for a linear 
array, the minimization problem with constraints can written: 



arg mm f{x s ,t* s ) = \ " x i? ~ (*J - ^f? 1 < * < N 

Propagation constraint : \x s — Xi\ — \t* s — t* \ 
Causality constraint : t* s < mini(t*) 



and Let suppose L 



so that X, — L is also a a solution of the 



minimization problem, ie. V f(X s — L) = 0) 
The Jacobian of / is written as: 



V/0r s ,t*) = 2 



E ( X s - Xi) ((X S - X. t f - (t* s - t*f\ 
i ^ ' 

Ein-tD^xs-xrf-itt-t*) 2 ) 



3 We say that the function is coercive 
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If X s being a critical point, this leads to two equations: 



Ei {*. - Xi) ((x s - x t f - (t* s - t*f) = (1) 

Ei (*: - ((^ - - (n - t*f) = o (2) 

Assuming that X s — L being also a critical point, this leads to two equations: 

' (x s - Xi - L) {[x a - x t - L) 2 - (t* s -t*-L) 2 )=0 (3) 
Ei (*: -t*s+L) ((x s -x t - Lf (t* s t* L) 2 ) = (4) 

By developing the equation (3) and by using the equation (1), then: 

(3) => Yl (*> - **) [(*- - - (*: - - 2L [(*« - **) - (*» - **)]] - x^-it: - t*y 

i i 

... +2L^\(x s -x % )-(t:-t:)] = Q 

i 

=> L (xs x l ?+L 2 J2 K*» *i) K - ( x s Xif-(t* s - t*) 2 

i i i 

... +Lj2(x s - Xi )(t*-t*)=0 



The set of constraints requires that the term (x s — Xi) z — (t* — t* ) 2 is null. 

i 

We get the simplified equation: 

LJ2 (x s - Xi) - (t* s - t*) =J2(x s - Xi) ((x s - Xi) - (t* s - t*)) 

i i 

In the cases where x s — Xi < for all i, the set of constraints is equivalent to 
(x s — Xi) — (t* s — t*) = 0. Thus, if one assumes that x s — Xi < for all i, i.e that 
the source is outside the array convex hull (a segment), we find that previous 
implications are equivalences and thus that equation (3) is verified. Operating 
in the same manner for the equation (4), we obtain the following equations: 



(4) -*:)[( 



(*: - t*f - 2L [(x s - Xi ) - (*: - *?)]] (x s - xi) 2 -(t* s - 1*) 



...-2L 2 J2\ixs-x t )-(t:-t*)]=0 



- 2L J2(t* - K) (x s - xt) - 2Lj2(n -t:f + Lj2(x s - x t f - - t*f 



7.1 Degeneration line for a linear antenna array 
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...-2L 2 Y,(xs-x t )-(t* s -t*)=Q 

i 

Using the set of constraints as above, we obtain the following equation: 

LJ2 (*a - Xi) - K - t*) = J2 (*J - O ((*. - Xi) - (t; - t*)) 

i i 

The same analysis as above gives us the condition that the source is out of 
the antennas convex hull. This degeneration is an important point because it 
determines the convergence of minimization algorithms. In this case the problem 
of the reconstruction is ill-posed. 

The generalization of the previous calculation to higher dimensions is more 
delicate, insofar as there are infinitely many directions in which the source can 
move. The idea now is to translate the source, from its position rt, simultane- 
ously in all directions rt — rt and with the same distances. We define the unit 



vector on the direction source- antenna. It will be noted: e, 



The 



translation spatial direction thus defined, is given by the vector L — = 



Ei 



r s . Considering the reduced tem- 
poral variables, the wave required delay to traverse the distance induced by the 



translation 



L 



L 
V: 



. We write 



. Let V the vector of coordinates V = (~L , 
the first order optimality condition for the vector of M 4 : X s 

Vf(X s -V)= (E/i (X, - V)) ■ M ■ (X s — V) — M ■ (X s - V) ■ X t ) 
By introducing the condition V/ (X s ) = which implies that: 

(J2 fi (Xs)) -M-X S -M-(J2 fi (Xs) X^ = 

we obtain: 



V/ (X s — V) ~N (V T ■ M ■ V) ■ M ■ 



X x 



i 

2M (J2 (X s - XifM ■V-jfy-2(22(Xa- X l ) T M ■ v) ■ M ■ (X s — V) 
-(j2fi(Xs))-M-V 
According to the imposed form of the vector 7", then: 

V T ■ M ■ V = (~t ||^ J M (~t j~t J T = 
It remains then the following expression: 

V/ (X s - V) =2M (J2 (Xs - Xif M -V -X,) - 2(J2(X S - X t f M ■ v) ■ M ■ (X s - V) 



■M-V 
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The resolution of this equation should lead to an analytical expression for the 
topology of critical points. We failed to develop it, but we can already see 
that the explicit development leads to cross terms that will make simplifications 
difficult. Therefore, we have tried again an intuitive approach based on the 
numerical simulations presented section 3.3. 



8. Appendix 3 

8.1. Convex hull for a linear antenna array 

Let us consider the sub-array of the 3 aligned upper antennas presented in 
Fig. |J. The fi gure flOl shows the physical principle of the reconstruction of the 



External offline 
source 




Ground 



Figure 10: Scheme of the reconstruction problem of spherical waves for a ID array of antennas. 
For this configuration, the convex hull is the segment shown in red. 



Three situations must be considered: 

• the source located inside the array; 

• the source located outside the array but on the detector axis; 

• the source located outside this main axis. The latter corresponds to the 
typical problems encountered with of the man-made emitters located on 
the ground. 

The first situation leads to 2 half- lines cutting each other in a single point: the 



solution is unique (Figure 111 and the localization problem is well-posed. The 
source is unique and inside the line segment linking the nearest antennas to the 
source. This segment correspond to the convex hull within this geometry. We 
can also note that only the two antennas flanking the source then play a role in 
its localization. The problem writes: 



8.1 Convex hull for a linear antenna array 
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arg min f(x s ,t*) = - ^2((x a - x % ) 2 - (t* - t*ff 

i=l 

Propagation constraint : \x s — Xi\ = \t* — t*\ 
Causality constraint : t* < mini(t*) 
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Figure 11: Phase space representation in the case of a linear array of three antennas (shown 
as green squares located at x\ = —200 m, X2 = m, X3 = 200 m). The source is located at 
x s = 60 m when the instant of the emission is taken as the time origin (t s = Os). Because the 
source is outside this sub-array, the constraints on the positions of antenna 1 and 2 lead to 
the same equation t a = 60 — x s (black line). Equation for the antenna 3 (blue line) leads to 
i s = —60 + x s . The causality conditions restrict the initial lines to two half-lines (red lines). 
The source location (black star) is at the intercept of the both half-lines. 



About the source on-axis, but outside the convex hull, the arrival times 
between the antennas, are no longer related to the source position, but to their 
locations. Whatever their positions, the time differences remain constant (for 
equally spaced antennas), ft becomes impossible to distinguish between two 
different shifted sources by any length. The only relevant information lies in 



the direction of propagation of the wave (see figure 12 1. This result appears by 
a degeneration of solutions because all points located on the half-line starting 
from the first tagged antenna are solutions of the problem which is ill-posed. 
The source is outside the convex hull of the antenna array. 

On the configuration where source located outside this antenna axis (problem 
in two dimensions), the solving starts with: 

N 

,*\2\2 



f(x s ,y s ,e s ) = \ Etc*, - x *) 3 + y* - (*: - ^*) 2 ) 2 



arg mm ..„_.. , , } 

«=1 



Propagation constraint : (x s — 2^) 2 + y 2 g — (t* — t* ) 2 
Causality constraint : t* < mirii(t*) 



8.1 Convex hull for a linear antenna array 
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Source outside convex hull 




Figure 13: Reconstruction problem of spherical waves for a ID array of antennas with the 
source outside the convex hull. The local minima are located at the intersection of the cones. 

The constraint set reduces the problem of characterization of critical points 
to the search of the half-cones intersections induced by each antenna, in the 
3 dimensional phase space (x, y , t) and which presents a great similarity of 



constraints with the light cone used in special relativity (Fig. 13 1. Intersection 



of the half-cones, two to two, induces multiple critical points which are local 
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